"""
# Cahn-Hilliard Equation - Spinodal Decomposition and Phase Separation

This module implements the Cahn-Hilliard equation, a fundamental model for 
describing spinodal decomposition, phase separation with conserved order parameter,
and coarsening dynamics. It's widely used in materials science, polymer physics,
and binary fluid systems.

## Mathematical Background

The Cahn-Hilliard equation is:
```
∂φ/∂t = ∇·[M∇(δF/δφ)] = ∇·[M∇(f'(φ) - ε²∇²φ)]
```

Where:
- φ: Conserved order parameter (concentration difference)
- M: Mobility function M(φ) 
- f(φ): Free energy density (double-well potential)
- ε: Interface width parameter

## Standard Form
```
∂φ/∂t = M∇²[φ³ - φ - ε²∇²φ]
```

## Key Features

- Mass conservation: ∫φ dx = constant
- Energy dissipative: dF/dt ≤ 0
- Fourth-order PDE (requires two boundary conditions)
- Coarsening dynamics: L(t) ~ t^(1/3)

## Applications

- Binary alloy phase separation
- Polymer blend demixing  
- Tumor growth modeling
- Image inpainting and denoising
- Microstructure evolution

## References

- Cahn, J. W. & Hilliard, J. E. (1958). J. Chem. Phys. 28, 258.
- Bray, A. J. (1994). Adv. Phys. 43, 357-459.
- Novick-Cohen, A. (2008). Handbook of Differential Equations, Vol. 4.
"""

using LinearAlgebra, FFTW, Plots, Statistics, Printf
using Random  # For reproducible random initial conditions

export CahnHilliardSystem, run_cahn_hilliard, analyze_spinodal_decomposition, measure_coarsening_dynamics

"""
    CahnHilliardSystem

Configuration for Cahn-Hilliard equation simulations.

# Fields
- `mobility_type::Symbol`: Type of mobility (:constant, :degenerate, :nondegenerate)
- `mobility_params::NamedTuple`: Parameters for mobility function
- `epsilon::Float64`: Interface width parameter ε
- `domain_size::Tuple`: Spatial domain dimensions
- `grid_size::Tuple`: Grid resolution
- `boundary_condition::Symbol`: Boundary conditions (:periodic, :no_flux)
"""
struct CahnHilliardSystem
    mobility_type::Symbol
    mobility_params::NamedTuple
    epsilon::Float64
    domain_size::Tuple{Float64, Float64}
    grid_size::Tuple{Int, Int}
    boundary_condition::Symbol
end

"""
    CahnHilliardSystem(; mobility_type=:constant, mobility=1.0, epsilon=0.1,
                       domain_size=(2.0, 2.0), grid_size=(128, 128), 
                       boundary_condition=:periodic, kwargs...)

Create Cahn-Hilliard system configuration.
"""
function CahnHilliardSystem(; mobility_type=:constant, mobility=1.0, epsilon=0.1,
                           domain_size=(2.0, 2.0), grid_size=(128, 128),
                           boundary_condition=:periodic, kwargs...)
    
    if mobility_type == :constant
        mobility_params = (M0=mobility,)
    elseif mobility_type == :degenerate
        mobility_params = (M0=mobility, gamma=1.0)
    elseif mobility_type == :nondegenerate  
        mobility_params = (M0=mobility, M1=0.1)
    else
        throw(ArgumentError("Unknown mobility type: $mobility_type"))
    end
    
    # Merge with user parameters
    mobility_params = merge(mobility_params, NamedTuple(kwargs))
    
    return CahnHilliardSystem(mobility_type, mobility_params, epsilon,
                             domain_size, grid_size, boundary_condition)
end

"""
    mobility_function(φ, system::CahnHilliardSystem)

Compute mobility function M(φ).
"""
function mobility_function(φ, system::CahnHilliardSystem)
    if system.mobility_type == :constant
        return system.mobility_params.M0 * ones(size(φ))
        
    elseif system.mobility_type == :degenerate
        # M(φ) = M0(1 - φ²)^γ (vanishes at φ = ±1)
        p = system.mobility_params
        return p.M0 * (1 .- φ.^2).^p.gamma
        
    elseif system.mobility_type == :nondegenerate
        # M(φ) = M0 + M1φ² (always positive)
        p = system.mobility_params
        return p.M0 .+ p.M1 * φ.^2
    end
end

"""
    chemical_potential(φ, ∇²φ, system::CahnHilliardSystem)

Compute chemical potential μ = δF/δφ = f'(φ) - ε²∇²φ.
"""
function chemical_potential(φ, ∇²φ, system::CahnHilliardSystem)
    # Standard double-well potential: f'(φ) = φ³ - φ
    f_prime = φ.^3 .- φ
    
    # Chemical potential
    μ = f_prime .- system.epsilon^2 * ∇²φ
    return μ
end

"""
    initialize_cahn_hilliard(system::CahnHilliardSystem, initial_type::Symbol=:spinodal)

Initialize the concentration field with appropriate initial conditions.

# Arguments
- `initial_type::Symbol`: Type of initial condition
  - `:spinodal`: Random noise for spinodal decomposition
  - `:nucleation`: Localized nuclei in metastable state
  - `:periodic`: Periodic perturbation
  - `:step_function`: Sharp interface
  - `:custom_noise`: Custom noise with specified statistics
"""
function initialize_cahn_hilliard(system::CahnHilliardSystem, initial_type::Symbol=:spinodal; 
                                 mean_concentration=0.0, noise_amplitude=0.1, seed=nothing)
    nx, ny = system.grid_size
    Lx, Ly = system.domain_size
    
    # Set random seed for reproducibility
    if seed !== nothing
        Random.seed!(seed)
    end
    
    # Create coordinate grids
    x = range(0, Lx, length=nx)
    y = range(0, Ly, length=ny)
    X = [xi for xi in x, yi in y]
    Y = [yi for xi in x, yi in y]
    
    if initial_type == :spinodal
        # Random noise around mean concentration (inside spinodal region)
        φ = mean_concentration .+ noise_amplitude * (2 * rand(nx, ny) .- 1)
        
    elseif initial_type == :nucleation
        # Small nuclei in metastable background
        φ = mean_concentration * ones(nx, ny)
        
        # Add nucleation sites
        n_nuclei = 8
        for _ in 1:n_nuclei
            center_x = Lx * rand()
            center_y = Ly * rand()
            radius = min(Lx, Ly) / 20
            distance = sqrt.((X .- center_x).^2 + (Y .- center_y).^2)
            φ .+= 0.5 * exp.(-(distance / radius).^2)
        end
        
    elseif initial_type == :periodic
        # Periodic perturbation
        kx = 2π / Lx
        ky = 2π / Ly
        φ = mean_concentration .+ noise_amplitude * (cos.(kx * X) .* cos.(ky * Y))
        
    elseif initial_type == :step_function
        # Step function interface
        φ = ifelse.(X .< Lx/2, -0.5, 0.5) .+ 0.01 * randn(nx, ny)
        
    elseif initial_type == :custom_noise
        # Gaussian random field with specified correlation length
        correlation_length = min(Lx, Ly) / 10
        φ = generate_correlated_noise(nx, ny, Lx, Ly, correlation_length)
        φ .*= noise_amplitude
        φ .+= mean_concentration
        
    else
        throw(ArgumentError("Unknown initial type: $initial_type"))
    end
    
    # Ensure mass conservation if specified
    if mean_concentration != mean(φ)
        φ .-= mean(φ)
        φ .+= mean_concentration
    end
    
    return φ, X, Y
end

"""
    generate_correlated_noise(nx, ny, Lx, Ly, correlation_length)

Generate spatially correlated random noise using FFT.
"""
function generate_correlated_noise(nx, ny, Lx, Ly, correlation_length)
    # Create Gaussian noise in Fourier space
    kx = 2π * fftfreq(nx, nx/Lx)
    ky = 2π * fftfreq(ny, ny/Ly)
    
    Kx = kx' .* ones(ny)
    Ky = ones(nx)' .* ky
    K = sqrt.(Kx.^2 + Ky.^2)
    
    # Correlation function in k-space (Gaussian)
    correlation_k = exp.(-(K * correlation_length).^2 / 2)
    
    # Generate random phases
    noise_k = (randn(nx, ny) + im * randn(nx, ny)) .* sqrt.(correlation_k)
    
    # Transform to real space
    noise = real.(ifft(noise_k))
    
    return noise
end

"""
    create_fourth_order_operator(system::CahnHilliardSystem)

Create operators for the fourth-order Cahn-Hilliard equation.
"""
function create_fourth_order_operator(system::CahnHilliardSystem)
    nx, ny = system.grid_size
    Lx, Ly = system.domain_size
    
    if system.boundary_condition == :periodic
        # FFT-based operators for periodic boundaries
        kx = 2π * fftfreq(nx, nx/Lx)
        ky = 2π * fftfreq(ny, ny/Ly)
        
        Kx = kx' .* ones(ny)
        Ky = ones(nx)' .* ky
        K2 = Kx.^2 + Ky.^2
        
        laplacian_k = -K2
        return laplacian_k, :spectral
    else
        # Finite difference for other boundary conditions
        dx, dy = Lx/nx, Ly/ny
        return (dx, dy), :finite_difference
    end
end

"""
    apply_cahn_hilliard_operator!(φ, μ, φ_new, laplacian_op, method, system::CahnHilliardSystem, dt)

Apply one time step of the Cahn-Hilliard operator.
"""
function apply_cahn_hilliard_operator!(φ, μ, φ_new, laplacian_op, method, system::CahnHilliardSystem, dt)
    if method == :spectral
        # Forward FFT
        φ_fft = fft(φ)
        
        # Compute Laplacian of φ in Fourier space
        ∇²φ_fft = laplacian_op .* φ_fft
        ∇²φ = real.(ifft(∇²φ_fft))
        
        # Compute chemical potential
        μ .= chemical_potential(φ, ∇²φ, system)
        
        # Apply mobility
        M = mobility_function(φ, system)
        flux = M .* μ
        
        # Compute divergence of flux
        flux_fft = fft(flux)
        div_flux_fft = laplacian_op .* flux_fft
        div_flux = real.(ifft(div_flux_fft))
        
        # Update concentration
        φ_new .= φ .+ dt * div_flux
        
    else
        # Finite difference method
        dx, dy = laplacian_op
        nx, ny = size(φ)
        
        # Compute Laplacian using finite differences
        ∇²φ = zeros(nx, ny)
        apply_laplacian_fd!(φ, ∇²φ, dx, dy, system.boundary_condition)
        
        # Compute chemical potential
        μ .= chemical_potential(φ, ∇²φ, system)
        
        # Apply mobility and compute divergence
        M = mobility_function(φ, system)
        flux_x, flux_y = compute_mobility_flux(M, μ, dx, dy)
        
        # Compute divergence of flux
        div_flux = zeros(nx, ny)
        compute_divergence!(flux_x, flux_y, div_flux, dx, dy, system.boundary_condition)
        
        # Update concentration
        φ_new .= φ .+ dt * div_flux
    end
end

"""
    apply_laplacian_fd!(u, ∇²u, dx, dy, bc_type)

Apply Laplacian using finite differences.
"""
function apply_laplacian_fd!(u, ∇²u, dx, dy, bc_type)
    nx, ny = size(u)
    
    for j in 1:ny, i in 1:nx
        if bc_type == :periodic
            i_prev = i == 1 ? nx : i-1
            i_next = i == nx ? 1 : i+1
            j_prev = j == 1 ? ny : j-1
            j_next = j == ny ? 1 : j+1
        else  # no-flux boundary conditions
            i_prev = max(1, i-1)
            i_next = min(nx, i+1)
            j_prev = max(1, j-1)
            j_next = min(ny, j+1)
        end
        
        ∇²u[i,j] = (u[i_prev,j] - 2*u[i,j] + u[i_next,j]) / dx^2 + 
                   (u[i,j_prev] - 2*u[i,j] + u[i,j_next]) / dy^2
    end
end

"""
    compute_mobility_flux(M, μ, dx, dy)

Compute mobility flux components.
"""
function compute_mobility_flux(M, μ, dx, dy)
    nx, ny = size(M)
    
    flux_x = zeros(nx, ny)
    flux_y = zeros(nx, ny)
    
    # Compute flux at cell faces (average mobility)
    for j in 1:ny, i in 1:nx-1
        M_avg = 0.5 * (M[i,j] + M[i+1,j])
        dμ_dx = (μ[i+1,j] - μ[i,j]) / dx
        flux_x[i,j] = M_avg * dμ_dx
    end
    
    for j in 1:ny-1, i in 1:nx
        M_avg = 0.5 * (M[i,j] + M[i,j+1])
        dμ_dy = (μ[i,j+1] - μ[i,j]) / dy
        flux_y[i,j] = M_avg * dμ_dy
    end
    
    return flux_x, flux_y
end

"""
    compute_divergence!(flux_x, flux_y, div_flux, dx, dy, bc_type)

Compute divergence of flux field.
"""
function compute_divergence!(flux_x, flux_y, div_flux, dx, dy, bc_type)
    nx, ny = size(div_flux)
    
    for j in 1:ny, i in 1:nx
        # x-direction divergence
        if bc_type == :periodic
            i_prev = i == 1 ? nx : i-1
            div_x = (flux_x[i,j] - flux_x[i_prev,j]) / dx
        else
            if i == 1
                div_x = flux_x[i,j] / dx
            elseif i == nx
                div_x = -flux_x[i-1,j] / dx
            else
                div_x = (flux_x[i,j] - flux_x[i-1,j]) / dx
            end
        end
        
        # y-direction divergence
        if bc_type == :periodic
            j_prev = j == 1 ? ny : j-1
            div_y = (flux_y[i,j] - flux_y[i,j_prev]) / dy
        else
            if j == 1
                div_y = flux_y[i,j] / dy
            elseif j == ny
                div_y = -flux_y[i,j-1] / dy
            else
                div_y = (flux_y[i,j] - flux_y[i,j-1]) / dy
            end
        end
        
        div_flux[i,j] = div_x + div_y
    end
end

"""
    calculate_cahn_hilliard_energy(φ, system::CahnHilliardSystem)

Calculate total free energy of the Cahn-Hilliard system.
"""
function calculate_cahn_hilliard_energy(φ, system::CahnHilliardSystem)
    nx, ny = size(φ)
    dx = system.domain_size[1] / nx
    dy = system.domain_size[2] / ny
    
    # Double-well potential: f(φ) = (φ² - 1)²/4
    potential_density = (φ.^2 .- 1).^2 / 4
    potential_energy = sum(potential_density) * dx * dy
    
    # Gradient energy: ε²/2 |∇φ|²
    if system.boundary_condition == :periodic
        # Use FFT for gradient calculation
        φ_fft = fft(φ)
        kx = 2π * fftfreq(nx, nx/system.domain_size[1])
        ky = 2π * fftfreq(ny, ny/system.domain_size[2])
        
        Kx = kx' .* ones(ny)
        Ky = ones(nx)' .* ky
        
        gradient_energy_k = 0.5 * system.epsilon^2 * (Kx.^2 + Ky.^2) .* abs.(φ_fft).^2
        gradient_energy = real(sum(gradient_energy_k)) * dx * dy / (nx * ny)
    else
        # Finite difference gradient
        ∇φ_x, ∇φ_y = finite_difference_gradient(φ, dx, dy)
        gradient_energy = 0.5 * system.epsilon^2 * sum(∇φ_x.^2 + ∇φ_y.^2) * dx * dy
    end
    
    total_energy = potential_energy + gradient_energy
    return total_energy, gradient_energy, potential_energy
end

"""
    finite_difference_gradient(φ, dx, dy)

Compute gradient using finite differences.
"""
function finite_difference_gradient(φ, dx, dy)
    nx, ny = size(φ)
    ∇φ_x = zeros(nx, ny)
    ∇φ_y = zeros(nx, ny)
    
    for j in 1:ny, i in 1:nx
        # Periodic boundaries
        i_next = i == nx ? 1 : i + 1
        i_prev = i == 1 ? nx : i - 1
        j_next = j == ny ? 1 : j + 1
        j_prev = j == 1 ? ny : j - 1
        
        ∇φ_x[i,j] = (φ[i_next,j] - φ[i_prev,j]) / (2*dx)
        ∇φ_y[i,j] = (φ[i,j_next] - φ[i,j_prev]) / (2*dy)
    end
    
    return ∇φ_x, ∇φ_y
end

"""
    run_cahn_hilliard(system::CahnHilliardSystem; t_final=1.0, dt=nothing,
                      initial_type=:spinodal, save_interval=0.1, verbose=true, kwargs...)

Run Cahn-Hilliard equation simulation.
"""
function run_cahn_hilliard(system::CahnHilliardSystem; t_final=1.0, dt=nothing,
                          initial_type=:spinodal, save_interval=0.1, verbose=true, kwargs...)
    
    if verbose
        println("🧪 Cahn-Hilliard Spinodal Decomposition")
        println("======================================")
        println("Mobility: $(system.mobility_type)")
        println("Interface width ε = $(system.epsilon)")
        println("Domain: $(system.domain_size)")
        println("Grid: $(system.grid_size)")
        println("Boundary: $(system.boundary_condition)")
        println()
    end
    
    # Initialize system
    φ, X, Y = initialize_cahn_hilliard(system, initial_type; kwargs...)
    nx, ny = system.grid_size
    
    # Check mass conservation
    initial_mass = sum(φ) * (system.domain_size[1]/nx) * (system.domain_size[2]/ny)
    if verbose
        println("Initial mass: $(round(initial_mass, digits=6))")
        println("Mean concentration: $(round(mean(φ), digits=4))")
    end
    
    # Create operators
    laplacian_op, method = create_fourth_order_operator(system)
    
    # Determine stable time step
    if dt === nothing
        if method == :spectral
            # Stability condition for fourth-order equation
            max_k4 = maximum(abs.(laplacian_op).^2)
            if system.mobility_type == :constant
                max_mobility = system.mobility_params.M0
            else
                max_mobility = maximum(mobility_function(φ, system))
            end
            dt = 0.1 / (max_mobility * max_k4 * system.epsilon^2)
        else
            dx, dy = laplacian_op
            min_dx = min(dx, dy)
            dt = min_dx^4 / (system.mobility_params.M0 * system.epsilon^2) * 0.1
        end
        dt = min(dt, 0.01)  # Upper bound
    end
    
    n_steps = Int(t_final / dt)
    n_save = max(1, Int(save_interval / dt))
    
    if verbose
        println("Time integration:")
        println("  dt = $(round(dt, digits=8))")
        println("  Steps: $n_steps") 
        println("  Save interval: $(save_interval)")
        println()
    end
    
    # Storage arrays
    times = Float64[]
    snapshots = Array{Float64,3}(undef, nx, ny, 0)
    energy_history = Float64[]
    mass_history = Float64[]
    characteristic_length_history = Float64[]
    
    # Work arrays
    φ_new = similar(φ)
    μ = similar(φ)
    
    # Save initial condition
    t = 0.0
    push!(times, t)
    snapshots = cat(snapshots, reshape(φ, nx, ny, 1), dims=3)
    
    # Calculate initial diagnostics
    total_energy, grad_energy, pot_energy = calculate_cahn_hilliard_energy(φ, system)
    push!(energy_history, total_energy)
    push!(mass_history, initial_mass)
    
    char_length = calculate_characteristic_length(φ, system)
    push!(characteristic_length_history, char_length)
    
    if verbose
        println("Initial diagnostics:")
        println("  Total energy: $(round(total_energy, digits=6))")
        println("  Characteristic length: $(round(char_length, digits=4))")
        print("Progress: ")
    end
    
    # Time evolution
    for step in 1:n_steps
        # Time step
        apply_cahn_hilliard_operator!(φ, μ, φ_new, laplacian_op, method, system, dt)
        φ, φ_new = φ_new, φ
        t += dt
        
        # Save data
        if step % n_save == 0
            push!(times, t)
            snapshots = cat(snapshots, reshape(φ, nx, ny, 1), dims=3)
            
            # Calculate diagnostics
            current_mass = sum(φ) * (system.domain_size[1]/nx) * (system.domain_size[2]/ny)
            push!(mass_history, current_mass)
            
            total_energy, grad_energy, pot_energy = calculate_cahn_hilliard_energy(φ, system)
            push!(energy_history, total_energy)
            
            char_length = calculate_characteristic_length(φ, system)
            push!(characteristic_length_history, char_length)
            
            if verbose
                progress = step / n_steps * 100
                print("$(Int(round(progress)))% ")
            end
        end
    end
    
    if verbose
        final_mass = mass_history[end]
        mass_error = abs(final_mass - initial_mass) / abs(initial_mass) * 100
        
        println("\n✅ Simulation completed!")
        println("Final mass: $(round(final_mass, digits=6))")
        println("Mass conservation error: $(round(mass_error, digits=4))%")
        println("Final energy: $(round(energy_history[end], digits=6))")
        println()
    end
    
    return (
        times = times,
        φ_final = φ,
        snapshots = snapshots,
        energy_history = energy_history,
        mass_history = mass_history,
        characteristic_length_history = characteristic_length_history,
        X = X,
        Y = Y,
        system = system,
        dt = dt
    )
end

"""
    calculate_characteristic_length(φ, system::CahnHilliardSystem)

Calculate characteristic length scale using structure factor.
"""
function calculate_characteristic_length(φ, system::CahnHilliardSystem)
    if system.boundary_condition != :periodic
        return 0.0  # Simple fallback for non-periodic
    end
    
    nx, ny = size(φ)
    
    # Remove mean for structure factor calculation
    φ_centered = φ .- mean(φ)
    
    # Calculate structure factor
    φ_fft = fft(φ_centered)
    S_k = abs.(φ_fft).^2 / (nx * ny)
    
    # Create radial average
    kx = 2π * fftfreq(nx, nx/system.domain_size[1])
    ky = 2π * fftfreq(ny, ny/system.domain_size[2])
    
    Kx = kx' .* ones(ny)
    Ky = ones(nx)' .* ky
    K = sqrt.(Kx.^2 + Ky.^2)
    
    # Find characteristic length (peak in structure factor)
    # Exclude k=0 mode
    K_flat = vec(K)[vec(K) .> 0]
    S_flat = vec(S_k)[vec(K) .> 0]
    
    if length(S_flat) > 0
        peak_idx = argmax(S_flat)
        k_peak = K_flat[peak_idx]
        characteristic_length = 2π / k_peak
    else
        characteristic_length = 0.0
    end
    
    return characteristic_length
end

"""
    analyze_spinodal_decomposition(results; verbose=true)

Analyze spinodal decomposition patterns and dynamics.
"""
function analyze_spinodal_decomposition(results; verbose=true)
    if verbose
        println("🔬 Spinodal Decomposition Analysis") 
        println("==================================")
    end
    
    φ_final = results.φ_final
    system = results.system
    times = results.times
    
    # Mass conservation check
    mass_history = results.mass_history
    initial_mass = mass_history[1]
    final_mass = mass_history[end]
    mass_error = abs(final_mass - initial_mass) / abs(initial_mass) * 100
    
    if verbose
        println("Conservation properties:")
        println("  Mass error: $(round(mass_error, digits=6))%")
    end
    
    # Energy dissipation
    energy_history = results.energy_history
    total_dissipated = energy_history[1] - energy_history[end]
    
    if verbose
        println("Energy dissipation:")
        println("  Total dissipated: $(round(total_dissipated, digits=6))")
        println("  Final energy: $(round(energy_history[end], digits=6))")
    end
    
    # Phase separation analysis
    φ_mean = mean(φ_final)
    φ_std = std(φ_final)
    
    # Calculate phase fractions
    positive_phase = count(φ_final .> 0) / length(φ_final)
    negative_phase = count(φ_final .< 0) / length(φ_final)
    
    if verbose
        println("Phase separation:")
        println("  Mean concentration: $(round(φ_mean, digits=4))")
        println("  Standard deviation: $(round(φ_std, digits=4))")
        println("  Positive phase fraction: $(round(positive_phase*100, digits=1))%")
        println("  Negative phase fraction: $(round(negative_phase*100, digits=1))%")
    end
    
    # Characteristic length evolution
    char_lengths = results.characteristic_length_history
    final_length = char_lengths[end]
    
    if length(char_lengths) > 2
        # Fit coarsening law L(t) ~ t^α
        valid_indices = char_lengths .> 0
        if sum(valid_indices) > 2
            log_times = log.(times[valid_indices][2:end])
            log_lengths = log.(char_lengths[valid_indices][2:end])
            
            if length(log_times) == length(log_lengths) && length(log_times) > 1
                # Linear regression for coarsening exponent
                A = [ones(length(log_times)) log_times]
                coeffs = A \ log_lengths
                coarsening_exponent = coeffs[2]
            else
                coarsening_exponent = NaN
            end
        else
            coarsening_exponent = NaN
        end
    else
        coarsening_exponent = NaN
    end
    
    if verbose
        println("Coarsening dynamics:")
        println("  Final characteristic length: $(round(final_length, digits=4))")
        if !isnan(coarsening_exponent)
            println("  Coarsening exponent: $(round(coarsening_exponent, digits=3)) (theory: 1/3)")
        else
            println("  Coarsening exponent: Not determined")
        end
        println()
    end
    
    return (
        conservation = (mass_error = mass_error,),
        energy = (dissipated = total_dissipated, final = energy_history[end]),
        phases = (mean = φ_mean, std = φ_std, 
                 positive_fraction = positive_phase, negative_fraction = negative_phase),
        coarsening = (final_length = final_length, exponent = coarsening_exponent),
        characteristic_lengths = char_lengths
    )
end

"""
    measure_coarsening_dynamics(results; plot_results=true)

Detailed analysis of coarsening dynamics and scaling laws.
"""
function measure_coarsening_dynamics(results; plot_results=true)
    println("📏 Coarsening Dynamics Analysis")
    println("===============================")
    
    times = results.times
    char_lengths = results.characteristic_length_history
    
    # Remove early time points and zeros for log-log analysis
    valid_mask = (char_lengths .> 0) .& (times .> 0.01)
    valid_times = times[valid_mask]
    valid_lengths = char_lengths[valid_mask]
    
    if length(valid_times) < 3
        println("Insufficient data for coarsening analysis")
        return nothing
    end
    
    # Fit power law: L(t) = L₀ * t^α
    log_t = log.(valid_times)
    log_L = log.(valid_lengths)
    
    # Linear regression
    A = [ones(length(log_t)) log_t]
    coeffs = A \ log_L
    log_L0, alpha = coeffs
    
    # Calculate R²
    L_fit = exp(log_L0) .* valid_times.^alpha
    ss_res = sum((valid_lengths .- L_fit).^2)
    ss_tot = sum((valid_lengths .- mean(valid_lengths)).^2)
    r_squared = 1 - ss_res/ss_tot
    
    println("Coarsening law: L(t) = $(round(exp(log_L0), digits=4)) * t^$(round(alpha, digits=3))")
    println("Theoretical exponent: 1/3 ≈ 0.333")
    println("R² = $(round(r_squared, digits=4))")
    
    # Compare with theoretical prediction
    theoretical_alpha = 1/3
    relative_error = abs(alpha - theoretical_alpha) / theoretical_alpha * 100
    println("Relative error: $(round(relative_error, digits=1))%")
    
    if plot_results
        # Create log-log plot
        p1 = plot(valid_times, valid_lengths, marker=:circle, label="Simulation",
                 xscale=:log10, yscale=:log10, xlabel="Time", ylabel="Characteristic Length",
                 title="Coarsening Dynamics", linewidth=2)
        
        plot!(valid_times, L_fit, label="Fit: t^$(round(alpha, digits=3))", 
              linestyle=:dash, linewidth=2)
        
        # Add theoretical line
        L_theory = valid_lengths[1] * (valid_times / valid_times[1]).^theoretical_alpha
        plot!(valid_times, L_theory, label="Theory: t^(1/3)", 
              linestyle=:dot, linewidth=2, color=:red)
        
        # Energy evolution
        p2 = plot(results.times, results.energy_history, xlabel="Time", ylabel="Energy",
                 title="Energy Evolution", linewidth=2, legend=false)
        
        coarsening_plot = plot(p1, p2, layout=(2,1), size=(600, 500))
        
        return coarsening_plot, (exponent=alpha, r_squared=r_squared)
    end
    
    return (exponent=alpha, r_squared=r_squared)
end

# Demonstration examples
if abspath(PROGRAM_FILE) == @__FILE__
    println("🧪 Cahn-Hilliard Equation - Research Examples")
    println("=" * 50)
    
    # Example 1: Classical spinodal decomposition
    println("\n1. Classical Spinodal Decomposition")
    system1 = CahnHilliardSystem(mobility_type=:constant, mobility=1.0, epsilon=0.05)
    results1 = run_cahn_hilliard(system1, t_final=2.0, initial_type=:spinodal,
                                mean_concentration=0.0, noise_amplitude=0.1)
    analysis1 = analyze_spinodal_decomposition(results1)
    coarsening_plot1, coarsening_data1 = measure_coarsening_dynamics(results1)
    
    # Example 2: Degenerate mobility
    println("\n2. Degenerate Mobility System")
    system2 = CahnHilliardSystem(mobility_type=:degenerate, mobility=1.0, 
                                 epsilon=0.03, gamma=2.0)
    results2 = run_cahn_hilliard(system2, t_final=3.0, initial_type=:nucleation)
    analysis2 = analyze_spinodal_decomposition(results2)
    
    # Example 3: Large-scale coarsening study
    println("\n3. Large-Scale Coarsening Study")
    system3 = CahnHilliardSystem(mobility_type=:constant, mobility=0.5, epsilon=0.02,
                                domain_size=(4.0, 4.0), grid_size=(256, 256))
    results3 = run_cahn_hilliard(system3, t_final=5.0, initial_type=:custom_noise,
                                correlation_length=0.2)
    analysis3 = analyze_spinodal_decomposition(results3)
    coarsening_plot3, coarsening_data3 = measure_coarsening_dynamics(results3)
    
    println("\n🎉 All Cahn-Hilliard examples completed!")
    println("Validated mass conservation and coarsening dynamics.")
    println("Observed coarsening exponents:")
    println("  Example 1: $(round(coarsening_data1.exponent, digits=3))")  
    println("  Example 3: $(round(coarsening_data3.exponent, digits=3))")
    println("  Theoretical: 0.333")
end